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Abstract. We focus on a particular connection between queueing and risk 
models in a multi-dimensional setting. We first consider the joint workload 
process in a queueing model with parallel queues and simultaneous arrivals at 
the queues. For the case that the service times are ordered (from largest in 
the first queue to smallest in the last queue) we obtain the Laplace-Sticltjcs 
transform of the joint stationary workload distribution. Using a multivariate 
duality argument between queueing and risk models, this also gives the Laplace 
transform of the survival probability of all books in a multivariate risk model 
with simultaneous claim arrivals and the same ordering between claim sizes. 

Other features of the paper include a stochastic decomposition result for 
the workload vector, and an outline how the two-dimensional risk model with 
a general two-dimensional claim size distribution (hence without ordering of 
claim sizes) is related to a known Riemann boundary value problem. 



1. Introduction 

There are several connections between queueing and risk models. A classical 
result is that the ruin probability in the Cramer-Lundberg risk model, in which the 
arrival process of claims is a compound Poisson process, is related to the workload 
(or waiting time) in an M/G/l queue with the same compound Poisson input. 
More precisely, denoting by (Rt)t>o the surplus process in the Cramer-Lundberg 
risk model, by r the time of ruin of this process and by (Vt)t>o the workload process 
in the corresponding M/G/l queue, one has P(r < t\Ro = u) = P(V t > u\Vq = 0); 
in particular, the probability of ruin ever occurring when starting at u equals the 
probability that the steady-state workload exceeds u. See, e.g., the nice geometric 
duality argument on page 46 of Asmussen and Albrecher [?] , or Rolski et al. [?] . 

However, also other ruin-related performance measures have a counterpart in 
queueing theory. By interpreting the interarrival times of the claims as service 
times of the corresponding queue and the claim sizes as interarrival times of the 
queue, the standard Cramer-Lundberg model is translated into a G/M/l queue. 
The time to ruin in the Cramer-Lundberg model is now related to a busy period of 
the corresponding queue, the deficit at ruin to an idle period and the surplus just 
before ruin to the sojourn time of the last customer in a busy period (see Frostig 
[?] and Lopker and Perry [?]). 
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In this paper our focus is on a connection between queueing and risk models in 
a multi-dimensional setting. In particular, we look at the joint workload process 
in a queueing model with parallel queues and simultaneous arrivals at the queues. 
Under the condition that, with probability 1, the service times of the customers 
arriving at the same time at the different queues are ordered (i.e., the customer in 
queue 1 has the largest service time, the customer in queue 2 the second largest 
service time, and so on) we are able to find the Laplace-Stieltjes transform of the 
joint stationary workload distribution in the different queues. Using a multivariate 
duality argument between queueing and risk models, this immediately gives the 
Laplace transform of the survival (non-ruin) probabilities in a multivariate risk 
model with simultaneous claim arrivals (and the same ordering property for the 
claim sizes of the simultaneous claims at the different books in the model) . 

Queueing models with parallel queues and simultaneous arrivals are also of- 
ten called fork-join queues. These models have many applications in computer-, 
communication- and production systems in which jobs are split among a number of 
different processors, communication channels or machines. Clearly, the queues in 
these models are dependent due to the simultaneous arrivals. In general this makes 
an exact analysis of the model very hard. Only in the case of two queues, exact re- 
sults are available (see, e.g., Flatto and Hahn [?], Wright [?], Baccelli [?], De Klein 
[?] and Cohen [?]). We will come back to some of these exact results in Section [5] 
of the paper. For the model with more than two servers no exact analytical results 
are available in the literature. In this case, bounds and approximations for several 
performance measures have been developed, see e.g. [?, ?, ?]. 

Multivariate risk models with simultaneous claim arrivals have several applica- 
tions in the area of ruin theory. One example is provided by reinsurance models 
in which, whenever a claim arrives, several insurance companies pay a part of the 
claim. Another example would be a large insurance company with multiple lines 
of business, where correlated claims arrive at the various business lines. Albeit in 
a different area of risk management, analysis of the dependence between the sto- 
chastic asset processes of several counter parties is also one of the most challenging 
aspects in the field of credit risk. Especially, in a two-dimensional setting one has to 
study the joint asset process of an obligor and a guarantor in credit default swaps. 

Avram, Palmowski and Pistorius [?, ?] have studied the joint ruin problem for the 
special case of two insurance companies that divide between them both claims and 
premia in some specific proportions. In particular, they derive the double Laplace 
transform with respect to the two initial reserves of the survival probabilities of the 
two companies. Proportional claims are a special case of our ordered claims, and we 
show in Section |4] that their survival result indeed is a special case of our Formula 
(JSJ). One of the key observations in [?, ?] is that, due to the fact that the companies 
divide the claims in some specific proportions, the two-dimensional ruin problem 
may be viewed as a one-dimensional crossing problem over a piecewise linear barrier. 
Badescu, Cheung and Rabehasaina [?] have extended the two-dimensional model 
of [?, ?] by allowing, next to the arrivals of claims for which the two insurers 
divide the claim in some specific proportions, also extra arrivals of claims which 
are fully paid by one of the insurers (e.g., insurer I). They show that under some 
conditions also in this model the previously mentioned reduction to a I-dimensional 
problem still holds. However, in [?] the authors do not consider the double Laplace 
transform with respect to the two initial reserves of the survival probabilities of the 
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two companies (their main focus is on the Laplace transform of the time until ruin 
of at least one insurer). 

The remainder of the paper is organized as follows: In Section [5] we present 
our model in detail and we provide the multivariate duality argument. This du- 
ality argument allows a translation between results for the queueing model and 
results for the multivariate risk model. Section [3] is dedicated to the analysis of 
the 2-dimensional queueing model with ordered service times. After introducing 
the assumptions, we derive the Laplace-Stieltjes transform of the joint stationary 
workloads in the two queues and present a decomposition theorem for the station- 
ary workload in the two queues. In Section [5] we extend the results of Section [3] 
to the A"-dimensional queueing model. Section |4] is dedicated to relations to other 
models. We present connections with tandem and priority queues, but also with 
a reinsurance problem with proportional claim sizes. In Section [6] we discuss the 
case of a general two-dimensional service time (or claim size) distribution. We in- 
dicate that the two-dimensional workload problem has been solved in the queueing 
literature. The solution is very complicated; our ordered service times case is a 
degenerate case, but a case which has the advantage of a much more explicit solu- 
tion which offers more probabilistic insight - and a case that can be generalized to 
higher dimensions. Finally, Section [7] outlines possible further research directions. 

Among the main contributions of our paper, we mention an explicit result for the 
transform of the joint workload (respectively, of the joint survival probability) and 
its extension to the A-dimensional model. In addition, we mention the workload 
decomposition result. It seems to be new in this setting, although similar results 
- under the assumption of independent inputs - were obtained for parallel queues 
(cf. [?]). From a more abstract perspective, another contribution of our paper is 
that it strengthens the links between queueing and risk models, pointing out that 
certain results and methods in the literature (and in the present paper) for queues 
with simultaneous arrivals are of immediate use in the risk setting, and vice versa. 

2. Multivariate Duality 

We consider a A-dimensional risk process in which claims arrive simultaneously 
in the K branches, according to a Poisson process with rate A. The claim sizes in the 
A books are independent, identically distributed random vectors (Bn\ '), 
n > 1. In the sequel we denote with (B^\ B^ K ') a random vector with the same 
distribution as (B^ , ...,b[ K '). 

For the n th arriving claim vector, denote by A n the time elapsed since the arrival 
of the previous claim vector, so that the A n are independent and have an identical 
exponential distribution with parameter A. 

(i) 

Let R t , i — 1, A be A risk reserve processes with initial capitals Ui, premium 
rates cW and the same arrival instants a n , n > 1. We have A n — a n — er rl _i and 
(To = (no delay). Then 

n(t) 

(1) R® = u i + £>«^ - Bf) + c« (t - a n[t) ) , 

i=i 

where n(t) is the number of arrivals before t. Let rW(iti) = inf 
be the times to ruin. 
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In connection with the ruin process, we consider K parallel M/G/l queues with 
simultaneous (coupled) arrivals and correlated service requirements. As in the ruin 
setting, A n are the interarrival times of customers in the K queues and the vector 
(B^, ...,B( K >) denotes the generic service requirements. The speed of server i is 
denoted by c^ l \ meaning that server i handles units of work per time unit, 
i=l,...,K. 

Furthermore we denote by pi := XE(B^) the load of queue i, i = 1, ...,K and 
we assume that pi < c^ % ' , to ensure that all queues can handle the offered traffic. 
These conditions imply positive safety loading in the ruin setting. 

From the queueing perspective, let {V t ,...,v} ) be the workload vector at 
time t in the system or, if we consider the n th arrival epoch, this is the workload 
(Vn X \ . . . , Vn ) seen by the customers of the n th batch arrival. Remark that 
Vn = c^Wn \ with Wn^ the waiting time of the n th arrival in queue i. Under the 
stability conditions above, the vectors {v} X \ V t ') and {Vn , Vn K ^) converge 
in distribution to the steady-state joint workload at arbitrary epochs and at arrival 
epochs, respectively. Due to the PASTA property these vectors are equal. Similarly, 
the vector {W n l \ . . . , W n K) ) converges in distribution to the steady state waiting 
time. We denote the Laplace-Stieltjes transform (LST) of the steady-state workload 
vector: 

For the multidimensional ruin process denned in ([1]), consider a dual workload 
process with vf? , i = 1, K the workload seen upon arrival by the N th customer 
in K initially empty queues with the time reverted arrival process (the arrival 
epochs are the same for all the systems): 

a *n = Cjv-n+ij {A^ = A N - n+ i), n = l,...,N; 

service time of customer n at queue i: Bn — B^_ n+1 , n = 1, N (time reverted 
service time) (cf. [?]). 

The following lemma shows that the well-known duality result (cf. [?], p. 46) 
between the Cramer-Lundberg model and the M/G/l queue can be extended to the 
multivariate risk model and the queueing model with simultaneous arrivals. Here 
the connection is between the various possibilities to be ruined (i.e we may have 
ruin in all books or precisely in one, at least in one, etc.) The results below are 
presented for the case K = 2, but can be directly extended to the general case. 

Lemma 1 . The following identities hold: 

(a) > ui A Vf } > u 2 ) - {tW(«i) < a N Ar^{u 2 ) < a N } 

(b) [v£ } < Ui A Vf > < u 2 ) = {T«(m) > „Ar( 2 )(« 2 ) > a N ) 

(c) {V^ > Ul AV^ < u 2 ) = {tW( Ui ) < a N AT^(u 2 ) > a N } 

(d) {V^ < Mi A ' > u 2 ) = {rW(ui) > o-n A r (2) (tt 2 ) < cr N } 
The above relations are path-wise identities. 

Proof. The following identities hold for the cylinder sets: 

{Vjp > Ui} = {r^(ui) < a N }. 
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This follows directly from Asmussen and Albrecher ([?], page 46) for the 1-dimensional 
problem, and is a special case of the duality in Siegmund [?] . 

If we intersect the above identities, we obtain (a). (&) follows by intersecting 
their complements, and (c) and (d) by subtracting (b) and (a) respectively, from 
the complements of the above cylinder sets. This concludes the proof. □ 

If we let N — > oo in (6) of Lemma [T] we obtain the infinite horizon joint survival 
probability 



(2) lim P(l4 1) < "i A < u 2 ) = P(r (1) (ui) = oo A t (2) (u 2 ) = oo). 

TV— voo 

Denote the righthand side by £(ui , U2) ■ This is the joint survival function, for initial 
capital (ui,U2). By PASTA, we can replace the steady state workload at arrival 
epochs with the steady state workload at arbitrary epochs in @. 



Let 



£.(a,t) := J e- SXl - tx ^(x 1 ,x 2 )dx 1 dx 2 



be the Laplace transform (LT) of the joint survival function. Via (J2|), this is also 
the LT of the c.d.f. of the joint workload in steady state. By a simple integration 
by parts, we have the following relation with the LST of the workload: 

(3) Us,t) = ~Ms,t). 

st 



3. The analysis of the two-dimensional problem 

In this section we derive the transform of the joint steady state workload process 
of the two-dimensional queueing model with simultaneous arrivals, as introduced in 
Section [2]. We also present a probabilistic interpretation of the quantities involved 
in the formula of the joint workload. The results are of immediate relevance for the 
corresponding insurance problem, via the duality outlined in the previous section. 

Before we start with the analysis, we make the following simplifying assumption. 
Assumption 1. All premium rates, respectively all service speeds, are 1, viz., 
C W = ... = C W=1. 

The following observation shows that this assumption is not restrictive. If we divide 
all terms in the righthand side of (JlJ by c^', we arrive at a new risk model with 
initial capital Ui/c^' and claim size B^/c^ and unit premium rates. Similarly, 
in the corresponding queueing model the service times at queue i are also divided 
by cW and the service speeds are equal to 1. This will not change the n th waiting 
time W$ at queue i, but the workload Vn^ at the n th arrival epoch is divided by 
. Also the times to ruin are preserved, hence the identities in Lemma [1] from 
the previous section remain unchanged. 

The LST of the joint service time/claim size vector is denoted by 

0(M) :=E(e~ sB(1> - tB{2> ). 
Our key assumption is the following: 

Assumption 2. P(B< 1 ) > B^) = 1. In view of the above discussion, in the case 
of speeds c (i) our assumption would be P(£W/c (1) > B^/c^) = 1. 
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Remark 1. This model allows for a dedicated Poisson arrival stream into queue 1. 
Merging this separate arrival process with the simultaneous arrival process at queue 
1, the distribution of B^ will have an atom in 0, which is the probability that a 
dedicated Poisson arrival happens instead of a simultaneous one ( see Badescu et al. 
[?] for a reinsurance model with both dedicated and simultaneous arrivals). 

We are interested in the joint stationary distribution of the amount of work in 
the two queues 

<KM) :=E(e- sVil> - tvi2> ). 
This can be obtained in the following way. Consider the amount of work in 
queue i just before the arrival of customer n. We have the following recursion for 
the random variables (Vn , Vn ),n = 1,2,... 

(V^% O - (max(K« + B« - A n , 0), max^j 2 ) + - A n , 0)). 
Or, for the LST 

M^t) = E(e-^- iV ^), n = l,2,..., 
this gives after straightforward calculations 

ip n+ l(s,t) = x _]_ t (<p(s, t)ip n (s, t) - (j>(s, A - s)i/j n (s, A - s)) 
+ (4>(s, A - s)i^ n (s, A - a) - 4>{\ 0)V„(A, 0)) 

(4) + 0(A,O)^„(A,O). 

Under the stability condition p\ < 1, ip(s,t) := limn-,.^ ijj n (s,t) exists and 
(i-i^)^(M) = (x^-*±h) A -a) 

(5) + (l-jij)^(A,0WA,0). 

If we let A denote a generic interarrival time, then due to the PASTA property, 

(6) 0(A, 0)^(A, 0) = P(V (1) + < A) = P{V W = 0) = 1- Pl . 

This is the probability that queue 1 is empty at an arbitrary time instant. 

On the regularity domains of ip(s, t) and <fi(s, t): We remark that, because of the 
dependence P(B^ > B^) = 1, we can rewrite the transform of the joint service 
times as: 

J>(s, t) = Ee-ffl' 1 '-^")-^™ = . fa s + tl 

and this function is always regular in IZe s > 0, lZe(s + t) > 0. If we consider 
(B^\B^) subject to > B^ a.s., cj)(s,t) may not be regular beyond this 

domain. More precisely, if B^ has a heavy-tailed distribution, this implies that 
£?W is also heavy tailed because of the dependence structure. In this case 4>(s,t) 
cannot be extended beyond IZe s > 0, IZe (s + 1) > 0. Similar considerations hold 
for ip(s,t) because we must also have P(V( X ) > V^) = 1. 

By Lemma [1] in the Appendix, Vs with IZe s > 0, there is a unique i(s), well 
defined and analytic in Re s > 0, such that A</>(s, t(s)) = X — s—t(s). Hence (s, t(s)) 

is a zero 

of ^ _ Mf^ i n p )j w hich is in the regularity domain of ip(s, t). Then 
the righthand side of ([5|) is also zero, i.e. 

(7) \t{s)(j){s, A - s)ip{s, A - s) = -s(A - t(s) - s)0(A, 0)^(A, 0). 
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If we substitute this in ([5]) and use ((6]), we obtain 

s t(s)-t 



(8) V(M) = (1-Pi) 



s + t-X(l-(p(s,t)) t(s) 



The interpretation of the Rouche zero t(s). Assume that a customer that 
starts a busy period BP( 2 ^ in queue 2 demands work x in queue 2 and work x + y 
in queue 1. During the service time of this customer in the second queue, there 
are Poisson(Ax) arriving customers, each one of these generating an i.i.d. busy 
sub-period with the same distribution as BP^ 2 ' in queue 2. So if we denote with U 
the extra work in the first queue, at the end of a busy period in the second queue, 
and with U*(s) its Laplace- Stieltjes transform, we have the identity: 



U*(s) 



r°° r°° °° (\rr\ k 

/ / e -syy2^§-e- Xx [U*(s)] k dP(BW-BW <y,B& <x). 
J x =o J y =o fe=Q k. 

The powers of U*(s) correspond to the extra work contributions at the end of the 
busy sub-periods started during the service time of the first customer in the busy 
period BP^. We can rewrite the above identity as: 

(9) U*(s) = 4>(s, A[l - U*(s)]) = <f>(s, A[l - U*(s)] - s). 
Comparing this with the equation in Lemma[T]in terms of s + t), we have: 

\4>(s,s + t{s)) = X - (s + t(s)) 
\j>(s,\[l-U*(s)]) =XU*{s). 

We may assume w.l.o.g. that P(B^ > B^) > 0, otherwise the two queues 
are a.s. identical, which is not interesting. Then it follows that the real part of 
A(l — U*(s)) is positive, and we must have s + t(s) = A(l — U*(s)) because the 
solution obtained in Lemma [T] is unique in the region Tie (s + t) > 0. We have thus 
proved: 

Proposition 1. The relation between t(s) and the transform of the extra workload 
in queue 1 at the end of a busy period in the shortest queue is 

(10) \U*(s) = A — (s + t(s)). 

The transform of the joint workload in the two systems becomes 

s + t- X(l-U*(s)) s 



V>M) = (i-pi) 



s + t- A(l -</){s,t)) s - A(l - U*{s)) 



The workload decomposition. Based on Proposition[TJ we show that the steady- 
state workload decomposes into an independent sum of a modified workload and an 
additional term, which represents the steady-state workload in a classical M/G/l 
queue. 

We start the joint workload process and let it run until the end of each busy 
period in the queue with the smallest workload. At this random time instant, we 
remove the extra content in queue 1, which has the largest workload of the two. Let 
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us denote this modified joint workload process as (V^, V^ 2 )). Then at the arrival 
instants of customers in the two queues, the recurrence relation holds: 

,>(i) t/ (2) v _ / + B n 1} - A n , K (2) + B n 2} - A n ), if A n < VP + B (2) 



\ v n+l> v n+l) " 



(0,0), if An > VP + B { n 2) . 



Remark that marginally, the shortest queue evolves unchanged. 

If we have ergodicity then in steady state, the above recurrence becomes: 

n/« T/(2h ± J + B{1) - A > y{2) + B(2) - A )i iiA < y{2) + B(2) 
(V ' K > ~ \ (0,0), if A>V^+B( 2 l 

Here and in the following, = denotes equality in distribution. If we rewrite 
this in terms of LST's, we obtain the following functional equation for tp(s,t) := 
Ee-'* m - tvW : 

(1 - ^ M . )i>(s, i) = (1 - 92) - , A M b, A - s)ct>{s, A - s), 
A — s — t A — s — t 

where 1 - p 2 = V(V^ = 0). 

Now follows a similar analysis as for ^(s, t). We already know from the Rouche 
problem that t(s) from Lemma [T] is a zero of (1 -*tLi&). We also have > F (2) 
a.s. (even if we take out the extra workload at the largest queue at the end of 
each busy period, is still at least as large as in the long run), therefore 
(s, i(s)) is in the regularity domain of ip(s,t) and therefore, at the point (s,t(s)), 
the right-hand side of the above identity is equal to zero: 

$(a, A - S )0( S , A - s) = (1 — p2) X ~ 3 ~ t{s) . 
Substituting back in the original identity, yields: 

This is a 2-dimensional Pollaczek-Khinchine type of representation. From an 
analytic point of view, the role of the numerator is to cancel the unique pole of the 
denominator in the region Tie (s + 1) > 0. 

Substitute t(s) from Proposition [1] and ip from (fTTj) into ([5]): 

(12) VM) = ^— ^ — * ^( a ,t). 

1 - Pi s - A[l - C/*(s)J 
We can now state the main result: 

Theorem 1 (Work decomposition). In steady state, we have the following repre- 
sentation of the joint workload at the two queues as an independent sum: 

(VW,VW)±(VW,VW) + (VW'\0), 

where yW' 1 is the workload in an independent, virtual M/G/l queue with arrival 
rate A and service requirements distributed as U , the extra workload at the end of 
a busy period BP*?' in the shortest queue. 
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Proof. It suffices to remark that the factor 

' _ =Ee ^ (1) ' 1 

l-p 2 s-X[l-U*(s)} 

in (fT2"j) is the Pollaczek-Khinchine formula for the transform of the workload in 
the virtual M/G/l queue with service time distribution U. This virtual queue 
is obtained by contracting the busy periods in the initial shortest queue, so that 
an arrival in the virtual queue happens at the end of this busy period and the 
interarrival time is then the idle period in the initial queue, and so is exponentially 
distributed. 

To see that indeed ^~ Pl is the atom of at 0, differentiate the identity for 

17* (s) in ©: 

ECU) = ~^(s, A(l - U*(s) - «)),,=„ = E(BW - B^) + XEB^ECU) 
as 

so that 1-XE(U) = □ 

4. Relation with other models 

In this section we point out how the results of the previous section are related 
to results for a risk model with proportional reinsurance, a particular tandem fluid 
model and with a particular priority queue. We start by showing that ([8]) generalizes 
a result obtained in [?], for the risk setting. 

The case of proportional reinsurance. In [?] the joint reserve process (R^\ R^) 
is of the form: R^(i) — m + c^'t/Si — S(t). Here S(t) is a common Compound 
Poisson input process with generic claim sizes a and are the premium rates. 
The claims are being divided in fixed proportions <5j, respectively. 

To bring this closer to our setting in Section [3j normalize the income rates: i.e. 
we consider -^R^) with pi = ^— . The assumption in [?] is that p\ > p 2 , 

which means that, in our notation, the claim sizes are := —a < —a =: B^ 2 \ 

Remark that the inequality between the B^'s is reversed here (which means the 
role of the arguments in our transforms is interchanged, especially the Rouche zero) . 
Let us recall the main formula in [?] (Formula (23)): 

/ -i o \ , / \ n 2 (0+y q+p- q + (q(pi - P2)) 

(13) ,_r<2) [p,q) = —, — -, — : — ; 7 tt y\ • 

p(K 1 {p + q) - q(pi -P2)) q-q + {q{pi -&)) 

The relation between the ruin times of (i? (1) , i? (2) ) and (— — i? (2) ) is 

TJ_ R (l) J_ R (2)(ui,U 2 ) = T RW M (2)(piUl,p2U 2 ). 

Pl ' P2 

Hence the relation to the LT coordinates used in ([3]) is s = p±p, t — p 2 q. From 
this, the relation between the LT of the survival functions becomes after a change 
of variables: 

( 14 ) i>*±RM ,-L-B,m{s,t) = -^—ip*R(i), R m(p,q)- 

pi 'P2 p x p 2 

• Hi (a) is the Laplace exponent of the Compound Poisson process with drift 
Pi per unit time. This means 

Ki(a) = piO, - A(l - Ee- aa ). 
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Because of the linear dependence between the B^'s, their LST has the 
form Ee- sS<1) - tB<2) = <j>(s, t) =: tj> BW (s + %t). 

• q + (q) is the largest root of the equation Ki(a) = q. Then q + (q(pi — P2)) 
solves: 

Pl a-\(l-Ee- a ^ B(1) ) = q( Pl ~p 2 ). 
Remark that if we set a = p + q, the above becomes: 
P\P + P2q - A(l - <j) B{ i) (pip + piq)) = 0, 

or, written in the (s, /^-coordinates, this becomes the equation in LemmaQ] 
(with s and t interchanged). Hence the relation between the zeroes in the 
two notations is: s(t) = pi(a — q) = Pi[q + {q(pi — P2)) — q]- 
The constant k(0+)' = P2 — XEB^ = p 2 (l — P2) is the probability that 
the queueing system is empty in steady state (now the second queue has a 
higher workload). 

In conclusion, (| 13[) written via (TT41 and © in the (s, t) coordinates becomes 
Formula ©: 

M t s) = S[l ~ P2) S ~ S{t) 

with the arguments s and t interchanged. 

Relation with work on tandem fluid queues. We now show that the workload 
model with ordered service times is equivalent with a particular tandem fluid queue. 
That is a model of two queues in series, in which the outflow from the first queue 
is a fluid, i.e., there is continuous outflow when the server is working (instead of 
customers leaving one by one). Such tandem fluid queues have been studied by 
various authors, see in particular [?]. Consider the following two-station tandem 
fluid network with independent compound Poisson input at the two stations (with 
arrival rate Ai and Laplace-Stieltjes transform of the service times B*(-),i = 1,2). 
Then Theorem 4.1 of Kella [?] gives the Laplace-Stieltjes transform of the steady- 
state fluid levels W\ and W2 in the two nodes: 

(15) iMax, a,) = E (e— — ■) = J 1 "*"^ ' ~ " ^ 

<j>i[a\) - <pi (572(0(2)) 0-2-772(02) 

with 

• Pl = X l E(B l ), 

• 4>i(ai) = ai - 771(011), 

. r) i (a i )=\ i (l-BZ(a i )), 

• 7)2 ("2) the solution of 0i (772(0:2)) = 772 (a 2 )- 

Alternatively, the last relation can also be formulated as: 7)2(02) is the solution of 

Ai-B 1 (f/ 2 (a 2 )) + X 2 B 2 (a 2 ) = Ai + A 2 - i) 2 (a 2 ). 

This system is related to our model with arrival rate A = Ai + A2 and Laplace- 
Stieltjes transform of service requirements 

^ s - f ) = irhr B ^ a + *) + rrr B ^ s )- 

Al + A2 Al + A2 



The corresponding notation is: B\ = B^ and B 2 = B^ — B^ 2 \ Here W± in the 
tandem model corresponds to the workload in the smallest queue in our model and 
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Figure 1 . Tandem fluid queue 



W\ + W% in the tandem model corresponds to the workload in the largest queue in 
our model. So we have 

il>{s,t) = E(e' sV ^)=E^e' s ^+ w ^- tw ^=i !w (s + t, S ) 

(1 - Pi - P2)s S + t- f) 2 {s) 

s + t-Xi(l-Bl(s + t))-X 2 (l-B^(s))' a -rials) 

Now remark that 

• The total traffic offered to the largest queue is pi+ p 2 , so indeed the factor 
1 — pi — pi in [?] corresponds to the factor 1 — pi in ©; 

. A(l - <t,{s, t)) = Ai(l - Bf(s + t)) + A 2 (l - B* 2 ( S )): 

• X(f>(s,t(s)) = XiBl(s + t(s)) + X 2 B* 2 (s) = Ai + A 2 - (s + t(s)), so indeed 
fj 2 (s) corresponds to our s + t(s). 

We conclude that flSJ coincides with Theorem 4.1 of [?] in the case of independent 
compound Poisson input. Kella's result is more general in the sense that he has 
Levy input instead of compound Poisson input. Our result is more general in the 
sense that we have dependent compound Poisson input. 

Relation with work on priority queues. As was already noticed in Kella [?], 
but also in several other places in the literature, the tandem fluid network described 
above is also related to a priority queue with preemptive resume priorities. Hence 
the same holds for our workload model. Consider the following model with two 
types of customers where customers of type-i arrive according to a Poisson process 
with rate A^ having service times with Laplace-Stieltjes transform B*(-),i = 1,2. 
Assume furthermore that customers of type-1 have preemptive resume priority over 
customers of type- 2. If we denote by Y\ and Y 2 the steady-state workloads in the 
two queues, then Yi and Y 2 are related to W\ and W 2 in the tandem fluid network. 
The Laplace-Stieltjes transform of the steady-state workloads in the two queues 
satisfies 

where again in our model we have to take arrival rate A = Ai + X 2 and Laplace- 
Stieltjes transform of service requirements 

Al + A 2 Ai + A 2 

We conclude that ([8]) also gives the Laplace-Stieltjes transform of a priority queue. 
Again our result is more general in the sense that we have dependent compound 
Poisson input (i.e., we can have arrivals of customers who have both low and high 
priority work) . 
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5. The ^-dimensional problem 

. In this section we consider the K-qnene system with simultaneous arrivals. We 
give the transform for the steady-state joint workload and we show that the de- 
composition in Theorem [T] extends to this case if we preserve the ordering between 
the service requirements/claim sizes. We use an iterative argument and for this 
purpose, the decomposition in Section [3] will be the starting point; the iteration 
step is essentially done with the help of Lemma [2] below as a work conservation 
identity. 

. We thus consider K parallel M/G/l queues, numbered 1 to K, respectively, with 
simultaneous (coupled) arrivals and correlated service requirements. The arrival 
process is again a Poisson process with rate A. The service requirements of suc- 
cessive customers at the K queues are independent, identically distributed random 
vectors (b£ } , . . . , B ( n K) ), n > 1. Let (#,...,#') be a generic random vector 
with the same distribution as (b[ 1} , b[ K) ). The LST of the service time/claim 
size vector is denoted by 

^s u ...,s K ):=E{e-^ BW -- s - s<f °). 

The essential assumption in the model extends Assumption 2 for the 2-dimensional 
problem: 

P(#>#>...>BW)=1. 

Furthermore we denote by pi := AE£?W, i = 1, . . . , K, the load of queue i and we 
assume that p± < 1 (hence pi < 1, Vi), to assure that all queues can handle the 
offered work. 

Remark 2. Like in the two-dimensional case (cf. Remark^), this model allows 
for a separate Poisson arrival stream into queue 1. Merging this separate arrival 
process with the simultaneous arrival process, the distribution of (B^ 2 \ . . . , Z?W) 
will have an atom in (0,...,0), which is the probability that a dedicated Poisson 
arrival happens instead of a simultaneous one. 

Similarly, the model allows for simultaneous arrivals at the first j queues only. 

This can be achieved by letting the distribution o/jB^ 1 ' BW) have an atom 

at (0, ...,0). 

The Laplace-Stieltjes transform of (V^, . . . , V^). Denote the Laplace-Stieltjes 
transform of the service time/claim size vector by 

^ Sl ,...,s K ):=E{e-^ BW --^ K) ). 

We have the ^-dimensional Lindley recursion for the random variables (V^ 1 "* , . . . , Vn K ^ ) : 

(V„% . . . , Vig) - (max(V« + - A n , 0), . . . , max(^W + - A n , 0)) n > 1. 

For the LST: 

M*i,---,«K)=E(e->rt 1) -"->* v ? e) ), n>l, 
the Lindley recursion gives after straightforward calculations: 



K 



ip n +i(si, ...,sk)=22 A „ E A j_ ~ [</> W (si> • • - , (si, • ■ • , Sj) 

3 = 1 

(16) - ^- 1 )( Sl ,..., Sj -_ 1 )^'- 1 )( Sl ,... )Si _ 1 )l +4> {0 H { ^ 
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where we used the following notation for simplicity: i()„ (s\ , . . . , Sk) '■— ipn(si, ■ ■ ■ , sk) 
and ipn^ := ipn(^, 0, . . . , 0), and 



ip^\si,...,Sj) := ip n (sx,...,Sj,X-^2si, O^^O ), for 1 < j < K - 1. 

1 ^ K—j—X arguments 

(fP> (si, . . . , Sj) is analogously defined for j = 0, . . . , K. By taking n. — > oo in (fT5|) . 
we obtain for ip(si,. . . , s K ) ■= lim^oo ip n (si,- ■ ■ , s K ), 

( 1 -^%rff)^.-'^) = E 

3=0 

(17) •0«(s 1 ,...,S,)^(s 1 ,...,S J ), 

with V W := hm ip^; and 0(°)t/>W = V(V^ + <A) = 1- Pl . 

Formula (I17I) has a simple recursive structure, and we can rewrite it as: 



(18) ^x-Vfa, sk-i) + (l - A ^-ff^' 0) ) . . . , 0). 

Denote by Cj := (l — A ^(si-_.sj.o,...,o) j ^(s^ . . . ) S j, 0, . . . , 0), and remark that 



?A(si, . . . , Sj ;, 0, . . . , 0) is the transform of the workload in the j-dimensional system 
obtained by ignoring the last (K — j) queues, j = 1, . . . , K. 

Proposition 2. The LST of the steady-state workload in the K > 3 systems is 

given by: 

(19) 

, _ (l-piOOSW-Sir) tt gj-gj 1 ~ Pi 5i 

E£i^-A(1-^i,..., S k)) M 1-PJ+i S j+ i l-p 2 S 2 > 
with Sj — Sj(si, Sj_i) i/ie unique solution of the equation 



\<f)(si,...,Sj,0,...,0) = A-^\ 



.; 

, s 
i=i 

wii/i IZe (si + ■ • • + Sj_i + Sj(si, . . . , Sj-i)) > 0, /or all j — 2, ... , K . 

Proof. The key remark is that sk is not among the arguments of the functions 
that appear in the righthand side of (jTT)) . 

From Lemma [T] applied to s = Si H h s^-i and t = sk, there exists a unique 

solution Sk = Sk(si, • ■ • , sk-i) of the equation 

K 

A0(si, . . . , s K ) = A - y^Sj; 

i=l 

such that Sk(si, . . . , s^-i) + X^i* s » nas positive real part. Hence the hyper- 
surface given by = Sk(s%, • ■ • , Sif-i) is contained in the regularity domain of 
ip(si, . . . , sk), and then the righthand side of (TT8l) must be zero. This gives the 
following relation for ^^"^(si, . . . , sk-i)- 
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AK-l) (a a \ i (K—l) I „ ^_ ^-^=1 s l )^Si, . . . ,SK-l,S K ) n 

<P (Si, . . . , Sk-1)W (Si, ■ ■ ■ ,sk-i) — ^ ^K-l- 

OK 

By substituting back into Equation (|18[) . we obtain the recursion 

n A - Sill 1 S i S K - S K 

with initial condition C2 = — (1 — pi) s 2 S2 , which follows from ([8]). From 

this, the formula in p^|) is obtained, after rearranging the factors. □ 

Interpretation of the Rouche root. It is worthwhile to change the coordinates: 
(si, s 2 , . . . , sr) (si, s 2 , . . . , SR-i,Yln=i s i)- We can rewrite 



(Si,...,SK 



) = Ee- s ^ B(1) - B<K ^--- SK -^ B(K ^- BlK ^-(^ S S< 



Let us denote it by (f>(si, . . . , sk-i, Ei=i s «)- This is the transform of the extra 
service time (relative to the shortest queue) in the first K—l queues, together with 
the shortest one. It turns out there is a connection between sk{si, ■ ■ ■ , sk-i) and 
the joint extra work in systems 1 to K — 1 at the end of a busy period in system K . 
Let us denote this extra work by (Z7i, U2, ■ ■ ■ , Uk-i), with LST U^(si, . . . , sk-i), 
and let F(xi, x 2 , . . . , xk) be the multivariate c.d.f. of 

- BW, . . . , B( K -V - B( K \BW). Then by a similar argument as the one 
leading to formula ©, U^(si, . . . , sk-i) satisfies the identity 

K-l 

U* K (s u sk-i) = fe~ £ ^^e- Xx « [U* K (si, s K -i)] n F(dxi . . . dx K ) 

J n=0 

(20) = 4>{si, sk-i, A[l - U* K (si, 8K-1)]). 
Comparing this with the identity for the Rouche root 

A — (si H h sk-i + S K ) = \<f>(si, . ■ • , sk-i,si H 1- sjc-i + Sk), 

gives the relation analogous to fj 10[) in Proposition 1 

(21) AC^(*i, . . . , = A - («i + • • • + s K _! + 
which follows because the Rouche root is unique. 

Let us fix our attention on the case K — 3 for the moment. Then identity (JT9J) 
becomes 

(22) 

(1 - p 3 )(S 3 - s 3 ) l-p 2 S 2 -s 2 l-pisi 



1p(si,S2,S 3 ) 



Si + S 2 + S3 - A[l - <f>(si,S2,S 3 )] 1~P3 S3 l-p 2 S 2 ' 



QUEUES AND RISK MODELS WITH SIMULTANEOUS ARRIVALS 



15 




V 31 



Figure 2. Work in the original system (left) and in the virtual 
system (right) 



Work conservation. We would like to give a probabilistic interpretation of (|2"2"1) . 

In order to achieve this, we start by considering the joint extra work in queues 1 
and 2 at the end of a busy period in queue 3. This has LST U£(si,S2) as input 
in a 2-dimensional system with simultaneous Poisson arrivals, which is obtained 
by contracting the busy cycles in queue 3. We call this the 2-dimensional virtual 
system. Remark that the inter-arrival times in the virtual system are precisely the 
idle periods in queue 3. 

For this construction, the key observation is that the steady-state extra work in 
the virtual queue 1 at the end of the busy period in the virtual queue 2 is the same 
as the extra work in the initial queue 1 at the end of the busy period in the original 
queue 2. In analytic form, let [/£ (si) be the LST of the extra work in the virtual 
system and J7| (si) be the LST of the extra work in the original system, see Figure 

m 

Lemma 2. 

u; (si) = u;( Sl ). 

Proof. We begin by remarking that the extra work (LA 1 )' 1 , U^' 1 ) in the first 2 
queues at the end of a busy period in queue 3 satisfies the a.s. inequality U^' 1 > 
[A 2 )' 1 . Since this is the input in the virtual system, from Proposition [IJ t/|(si) 
satisfies the identity ([9]) with C7|(si,S2) instead of </>(si,S2): 

(23) £/ 3 *(si,A[l - U2(si)] -ai) = £/ 2 *(si). 



At the same time, via (|20|) . U^(si,S2) satisfies 
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0(si, S2, A(l - [/* (si, s 2 )) - si - S2) = f7| (si, s 2 ). 
If we substitute this fixed point identity in ()23[) above, we have 



^( Sl ,A(l-f/ 2 *(s 1 ))-s 1 ,0) = C/ 2 *(s 1 ). 
On the other hand, this is also the identity ([9]) satisfied by U 2 (si), in the 2- 
dimensional system obtained by ignoring the last queue. Hence, from the uniqueness 
result in LemmalU U 2 is\) = U 2 (si) (See Figure[2]). This completes the proof. □ 

We can rewrite (19) using (21): 

si + s 2 + s 3 - A(l - U^(s 1 ,s 2 )) 



ip(st, s 2 , s 3 ) = (1 - p 3 ) 
(24) 



si + s 2 + s 3 - A(l - (t>(si,s 2 , s 3 )) 
l-p 2 s 1 + s 2 -\(l-U%(s 1 )) 1-pi si 



1 - p 3 si + s 2 - A(l - U*( Sll s 2 )) 1 - P2 Sl - A(l - U*( Sl )) ' 
Remark that the atom \Z P p l above is the conditional probability that queue 1 is 
empty, given that queue 2 is empty; and similarly for \Z_ P p 3 ■ m addition, the last 
factor in (|2"4")l is the Pollaczek-Khinchine representation for an M/G/l queue with 
service times having LST U 2 (s\). Now we are ready to give the main result of this 
section. 

Theorem 2. In steady state, the joint workload distribution decomposes as an 
independent sum: 

(VW, V^) £ (V-W.1, yWA F (3)) + (y(i),2 5 F (2),2 )0) + (F (i), 3j 0j 0) 

The first term in the sum represents the steady-state distribution of the modified 
joint workload process obtained by removing the extra work in the first two queues 
at the end of a busy period in the third queue. The second term is the workload 
in the first two queues obtained by removing the extra work in the first queue at 
the end of a busy cycle in the second queue. Finally the third term represents the 
workload in the virtual M/G/l queue with input distributed as the extra work in 
queue 1, at the end of a busy period in queue 2. 

Proof. Consider the modified work process that evolves in steady state as 

^(1)4 yW,l iV W) £ frKl),l + B (l) _ ^ y(2),l + B {2) _ A y(3) + B (3) _ A ^ 

if A < and V^' 1 , V&) = (0,0,0), else. 

By similar computations as the ones leading to Formula (jlip . we obtain 

7 1 \ /, x «i + s 2 + «3 - A(l - C/ 3 *(si,s 2 )) 

i>{si, s 2 , s 3 ) = (1 - p 3 )- 



si + s 2 + s 3 - A(l - 4>(si, s 2 ,s 3 )) ' 
This is the first factor in (l24l . For the second one, consider the following modified 
virtual workload process that evolves in steady state as 

( y(i),2 y(2),2 0) £ { (> (1) » 2 + A, + UW' 1 -A,0), if A < - 

\ (0,0,0), HA>V^ 2 - 
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with (LA 1 )* 1 , C/t 2 )' 1 ) the extra work vector in the first 2 queues at the end of a busy 
period in queue 3. Here we remove the excess workload in the virtual queue f at 
the end of the busy period in the virtual queue 2, which by Lemma [5] is the same 
as in the original system. In terms of LST's , this becomes 

7 , v I -Pi Sl+*2-A(l-E/J(*i)) 



1 - Pi si + s 2 - A(l - U*(s 1 ,s 2 ))' 

Finally, the third factor in (IM1) is the Pollaczek-Khinchine representation of the 
steady-state workload in the M/G/l queue with service time distributed as the 
extra work in queue 1 at the end of a busy period in queue 2. This ends the 
proof. □ 

These considerations can be iterated now for the general ii"-dimensional system. 

Corollary 1. The steady-state joint workload in the K systems decomposes into 
the independent sum 

{V [1 \...,V [K) ) = (V {1)tl ,...,V {K ~ l)>1 ,V {K) ) + (V (lh2 7 . . . ,V {K ^ 2) ' 2 ,V {K ^ 1) ' 2 ,0) 

+ ...+ (yW**- 1 , VW> K -\Q, . . . , 0) + (V^' K , 0, . . . , 0), 
where the jth term in the sum satisfies the identity in distribution (j—2,. . . ,K): 

(V-(i)j s y(2) J 5 . . . ; yiK-j) J > V {K-j+i) ,3 1 0j . . . > o) I (yWJ + U^'J- 1 - A, 

y{2),j + _ ^ f yiK-j+DJ _ B {K-3+l) _ A, 0, . . . , 0) , 

if A<V {K -i +1) 'i -B^-J+V, 

and (0, ... ,0) else. is the extra workload in queue i at the end of a busy 

period in queue (K — j + 1), for i > K — j + I. 



6. The general two-dimensional workload/reinsurance problem 

In this section we consider the general two-dimensional workload problem: pairs 
of customers arrive simultaneously at two parallel queues Q\ and Qi according to 
a Poisson(A) process, the nth pair requiring service times (B^\B^) with LST 
(f>(s,t). We are interested in the steady-state workload vector (V^yV^) with 
LST ip(s,t). By the duality that is exposed in Section^ ip(s,t) also is the Laplace 
transform (w.r.t. u\ and u 2 ) of the probability that both portfolios of an insurance 
company with simultaneous claims (Bn \ B^), with initial capital U\ and u 2 , will 
survive. 

In Section[3]we have determined ij}{s, t) for the special case that P(i?( 1 ^ > B 1 - 2 ^) = 
1. We now show how the general case - Bn^ and B^ having an arbitrary joint 
distribution - has been solved in the literature (with the solution of that special case 
emerging as a degenerate solution). We shall successively discuss the contributions 
of Baccelli [?], De Klein [?] and Cohen [?], who have treated the two-dimensional 
workload problem with simultaneous arrivals in increasing generality. Starting point 
in all those three studies is the following functional equation for tp(s,t), which 
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is derived by studying the 2-dimensional Markovian workload process during an 
infinitesimal amount of time At: 

(25) K(s, t)ip(s, t) = #i(s) + sip 2 (t), He s, t > 0. 
Here the so-called kernel K(s,t) is given by: 

(26) K(s,t) := s + t- X(l-cj)(s,t)), 
and 

(27) Ms) ■= ne- sVl (V2 = 0)], Mt) ■= E[e- tVa (7i = 0)], 
with (•) denoting an indicator function. 

Remark 3. In the special case of Section^ with ¥(B^ > B^) = 1, one has 
(t) = P(Vi = 0), because Vi cannot be positive when V\ — 0. It then remains to 
find ipi(s) . This is done by observing (cf. the appendix) that, for all s with Tie s > 0, 
there is a unique zero t(s) of the kernel, with He t(s) > He (— s). This immediately 
yields that ipi(s) — ~ tHS^O^ = ®)> w hich is readily seen to be in agreement with 

W- 

Equation {5p, which was obtained by studying the workloads at arrival epochs 
(i.e., the waiting times; by PASTA they have the same distribution as the steady- 
state workloads), looks slightly different from 125\) . but using it is readily seen 
that they are equivalent. 

Globally speaking, the essential steps in [?, ?, ?] are the following. 
Step 1: find a suitable set of zeroes (s,t), with Tie s > 0, He t > 0, of the kernel 
K(s,t), i.e., K(s,i) = 0. Because i^(s,t) is regular for all (s,t) with Tie s,t > 0, 
one must have for all these zeroes: 

(28) t>i(s) = -#2®. 

It is further observed that tpi(s) is regular for Tie s > 0, continuous for He s > 0, 
and that ip2(t) is regular for He t > 0, continuous for He t > 0. 
Step 2: formulate a boundary value problem for ipi(s) and ip2{t). There are various 
types of boundary value problems, like the Riemann and the Wiener-Hopf bound- 
ary value problems. Typically, they ask to determine two functions Pi(-) and PaO)) 
which satisfy a relation on a particular boundary B, while Pi(-) is regular in the 
interior B + and P2O) is regular in the exterior B~ . B could be the unit circle 
(Riemann boundary value problem), or the imaginary axis (Wiener-Hopf bound- 
ary value problem; B + now is the left-half plane). We refer to Gakhov [?] and 
Mushkelishvili [?] for excellent expositions of such boundary value problems and 
their variants, like the boundary value problem with a shift. The latter occurs in 
the approach of De Klein [?], see below. 

Step 3: solve the boundary value problem for ^i(-) and ^(O with boundary B. If 
B is a smooth closed contour that is not a circle, the use of a conformal mapping 
from B to the unit circle C is required to arrive at a Riemann boundary value 
problem for the unit circle, the solution of which can be found in [?, ?]. Thus one 
obtains ibi(s) and ip2{t) inside certain regions; subsequently, one may use analytic 
continuation to find them in He s, t > 0. Finally, ip(s, t) follows from ([25]) . 

Remark 4. Application of the boundary value method in queueing theory was pi- 
oneered by Fayolle and Iasnogorodski in [?]. They used this method to analyze 
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the joint queue length process in two coupled processors, viz., two M/M/l queues 
which operate at unit speeds when the other queue is not empty, but at different 
speeds when the other queue is empty. The method was subsequently developed in 
[?] for a large class of two-dimensional random walks; various queueing applications 
were also discussed in [?] . See [?] for a survey of the method in queueing theory, 
and see [?, ?] for two monographs which have further developed the theory of two- 
dimensional random walks. Part IV of [?] explores the analysis of N -dimensional 
random walks with N > 2. Results for N > 2 are very limited, and it seems fair to 
conclude that the boundary value method is, apart from a few special cases, restricted 
to two-dimensional random walks. 

Remark 5. We strongly believe that the boundary value method also has a large 
potential in the analysis of two-dimensional risk models. Due to the duality between 
the reinsurance model and the 2-queue model with simultaneous arrivals, the pub- 
lications [?, ?, ?] are of immediate relevance to the reinsurance problem. These 
publications seem unknown in the insurance community (see, e.g., Chan et al. [?], 
who pose the two-dimensional risk problem and stop at Equation &25]) (where [?, ?, ?] 
begin). They have remained largely unnoticed even in the queueing community, per- 
haps because of their complexity and because [?] and [?] did not appear in the open 
literature. For these reasons, we now successively expose the approaches in [?], [?] 
and [?] at some length. 

The approach of Baccelli [?] 
Baccelli [?] restricts himself to the case of exchangeable (B^\ B^), i.e., P(B^ < 
x,B {2) < y) = P(B (1) < y,B {2) < x), or equivalently, <j>(s,t) = <f>(t,s). We briefly 
review the three steps mentioned above. 

Step 1 in [?] is as follows. Consider zero pairs (s, t) = (g + iu,g — iu) of kernel 
K(s, t), with u 6 1Z and with g — g(u) the unique zero in Tie g > of 

2g = A(l - (f>(g + iu,g~ iu)). 

Using the exchangeability, it can be shown that this unique zero is real and non- 
negative, while g(— u) — g(u), u GlZ. 

Step 2. Consider the arc A = {s : s = g(u) + iu, u £ R}, with g{u) the zero 
defined above. This is a smooth arc, located in the right half-plane. Baccelli finds 
a conformal mapping p(-) of the interior C + of the unit circle C onto A + , the 
'interior' of A located on the right of A, and a conformal mapping q(-) of C~, the 
exterior of the unit circle, onto A + ; their limits on C are denoted by p + (z) and 
q~(z), which are each other's complex conjugates because of the exchangeability. 
Noticing that p + {— 1) = q~{— 1) = 0, he multiplies both sides of (|28p with 1 + z. 
This yields (divide both sides of ([25} by si): 

(29) (l + z) +l s. = -(! + *) —, -s. , \z\ = 1- 

p+(z) q (z) 

Because of the regularity properties of the conformal mappings and of V , i( s ) an d 
ip2(t), IZe s, t > 0, one now arrives at a simple boundary value problem: we have 
([29| for \z\ = 1, while the left-hand side of (|29|) is regular for \z\ < 1, and the 
right-hand side is regular for \z\ > 1. 

Step 3. The solution of this problem immediately follows from Liouville's theorem, 
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Cf. [?] p. 85: 

Baccelli [?] shows that 7 = —(5, and determines the remaining unknown constant 8 
by normalization. Having thus determined ipi( s ) f° r s € ne uses analytic con- 
tinuation to obtain ifti(s) in the whole right half-plane; similarly for ip2 (t)- Finally, 
substitution in (|25[) determines i/)(s,t). 

The approach of De Klein [?] 

De Klein [?], pp. 119-168, studies the general case of an arbitrary joint distribution 
of flW and B {2 \ 

Step 1 in [?] is as follows. He considers the same zero pairs as Baccelli (also 
suggesting another set of zero pairs on p. 132). g(u) is no longer necessarily real, 
but for all real u there still is a unique zero g{u). 

Step 2. De Klein subsequently considers the simple, smooth arcs A\ = {s : s — 
g(u)+iu, u £ R} and A2 — {t : t = g(u)—iu, u £ R} in the right half-plane. Notice 
that A\ and A2 are each other's complex conjugates in the exchangeable case of 
Baccelli, but not in De Klein's more general case. De Klein now uses the (unique) 
one-to-one mapping t = u>2(s) from A\ onto A2 (with inverse ui(t)) determined by 
the fact that, Vs £ A\, (s,Ld2(s)) is a zero pair of the kernel. Similarly, V< £ A2, 
(uii (t),t) is a zero pair. Hence the following must hold: 



(30) ^( Wl (t))=-^ 2 (t),t6A 2 . 

In addition, one has the regularity properties of the functions ipi(-) and "02(0 which 
were listed below (|2"8l . Determination of functions ipi(-) and - 2 (O with these 
regularity properties and satisfying (|30[) is a so-called shift problem, a boundary 
value problem with a shift (cf. Sections 17 and 18 of [?]). 

Step 3. Gakhov [?] mentions two methods to solve such problems: (i) reduce the 
problem to a Fredholm integral equation of the second kind, and (ii) reduce the 
problem to an ordinary Riemann boundary value problem, by means of conformal 
mappings. De Klein [?] explores the first method in Section II. 4. 2 and the second 
in Section II. 4. 3. We concentrate on the first method. De Klein first translates the 
shift problem to one on a finite smooth closed contour, via the conformal mapping 
Q(z) = (with inverse z(£) = j^) that maps A. k onto smooth closed contours T,-, 
2 = 1,2; he then applies Gakhov's first method. He obtains the following Fredholm 
integral equation of the second kind for an unknown function G\{-) - which up to 
a constant equals logl^i (;?(•))}: 
(31) 

Gi(pi) = ^-.\ Gi(vi)[— V '^ Vl) , — ]d«i+JZi(pi), pi £ Ti, 

with Hi(-) some known function, c\ some point in the interior of Ti, and 1^2 («i) 
= ^(^(^^i))), v\ £ Ti. After having solved the integral equation (which can be 
done numerically in an efficient way, as shown by De Klein), one obtains V , i( s ) f° r 
s £ Ti, and then 02 (i) for t € T 2 via (|28p . The regularity of ipi(s) in the interior 
Ti subsequently allows one to obtain - 0i(s), s £ T x + , as a Cauchy integral; similarly 
for tp2(t), t £ T 2 + . By analytic continuation, Vi( s ) an d ^(t) are then also uniquely 
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determined in TZe s > and IZe t > 0, respectively. Finally, ijj(s,t) again follows 
from (p5|), 

De Klein also explores Gakhov's second method to treat the shift problem. How- 
ever, this reduction to a Riemann boundary value problem requires a conformal 
mapping that itself must be determined by solving another Fredholm integral equa- 
tion of the second kind. In Chapter II. 6 he extensively investigates the numerical 
solution of both integral equations by means of the Nystrom or quadrature method. 
He obtains, a.o., accurate results for the mean sojourn time of a customer pair, viz., 
the time until both customers of a pair have left the system. 

The approach of Cohen [?] 

Cohen [?], Part III, considers a very general class of two-dimensional workload 
processes. Basically, he combines the model with simultaneous arrivals and the 
coupled processors model. The two servers have speeds r\ and r-i if they are both 
non-idle, and speeds r^ and A 2 ^ when the other server is idle. Furthermore, he 
also allows the possibility of different joint service requirement distributions if a 
customer pair arrives when at least one of the servers is idle. Finally, he explicitly 
allows single arrivals next to simultaneous arrivals (cf. also [?]). Much of Part III 
of [?] is devoted to a detailed study of the ergodicity conditions and of the so-called 
hitting point process and hitting point identity of the workload process, hitting 
point referring to the first entrance point of one of the axes. 

In Chapter III. 4 he determines the steady-state joint workload distribution for 
a variety of cases. For us, the most relevant cases are treated in Sections III. 4. 9 
and III. 4. 10. Section III. 4. 9 treats the model of De Klein [?]. The same zero 
pairs are used (Step 1), and the same smooth closed contours T\ and T2; Cohen 
subsequently uses Gakhov's second method to arrive at a Riemann boundary value 
problem (Step 2). That boundary value problem actually is so simple that it can 
be solved straightforwardly by applying Liouville's theorem, cf. Baccelli's method 
above (Step 3); however, a conformal mapping is required, which is obtained as the 
solution of another Fredholm integral equation of the second kind. A nice feature 
in Section III. 4. 9 is that fais) an d fait), after normalization, are expressed as 
LST's of waiting time or workload distributions of special M/G/l queues (which 
are related to hitting points). 

Section III. 4. 10 treats the model of De Klein with the additional feature that 
there is coupling of the servers, of a rather special form: + -^§y = 1. This does 
not change the kernel K(s, t) (which only refers to the interior of the state space, 
with both servers active), so the same zero pairs and contours can still be used. 
However, it does change the right-hand side of ([^5)1 . and hence a slightly different 
Riemann boundary value problem must be solved. 

Remark 6. It should be observed that Baccelli [?], De Klein [?] and Cohen [?] all 
also solve the more complicated transient problem, of determining the joint time- 
dependent distribution of the two workloads. 

7. Conclusions and future work 

We have studied a multivariate queueing system, which is shown to correspond 
to a dual risk process with multiple lines of insurance that receive coupled claims. 
We find the LST of the multivariate workload distribution in the case in which 
the service requirements are ordered with probability one. Duality then yields the 
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Laplace transform of the survival probabilities. For general service requirement 
(resp. claim size) vectors the workload (resp. ruin) problem can be solved in 
the two-dimensional case, by solving a Riemann boundary value problem. For 
dimension K > 2, the problem seems analytically intractable in its full generality. 
That raises the need for approximations and asymptotics. It would in particular be 
interesting to obtain explicit multi-dimensional tail asymptotics of workloads and 
ruin probabilities, both for light-tailed and heavy-tailed service requirements (or 
claim sizes). Even for K = 2 queues, this is already quite challenging. Moreover, a 
wide range of different cases must be studied, giving rise to quite different techniques 
and results. Therefore we intend to devote a separate study to tail asymptotics. 

8. Appendix 

Lemma 1 (Rouche zero). For every s with IZe s > there exists a unique t = t(s) 
with IZe t(s) > TZe (— s), that satisfies the identity 

\(j){s,t) = A- (s + t). 

Moreover the function: s — > t(s) is analytic in TZe s > 0. 

Proof. For fixed s with TZe s > 0, let f(s + t) := A — (s + 1). Consider in the right 
half-plane the contour C made up from the semicircle with center at — s and radius 
R > 2A together with the line segment / := {— s + iw\w e [— R, R]}. We show that 
on this contour \\<j)(s,t)\ < \f(s + t)\. We can bound \(j)(s,t)\ by 

X\cP(s,t)\ = \\4>(s,S + t)\ < XE e -Res(B^-B^)-Re(t+s)B^ <x 

This holds everywhere in the domain of (j)(s,t) if — has positive mass on 
(0,oo). 

Now we bound \f(s + t)\. When (s + t) is on the semicircle (i.e \s + t\ = R > 2A), 
apply the triangle inequality to the triangle with vertices at 0, A, s + t, to find 
A — s — t\ > A. When (s + t) e /, by a similar argument we obtain |A — s — t\ > A, 
with equality only when s + 1 = 0. Hence on the contour C, \f(s + t)\ > A. We can 
now use Rouche's theorem to conclude that the equation A</>(s, t) = A — (s + 1) has 
a unique solution t(s) inside C, because the polynomial f(s + t) has only one zero 
inside C, at A. Letting R — V oo, proves the assertion. □ 
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